'''
Created on Nov 27, 2010

@author: oabalbin
'''
import subprocess
import snps.samtools_commands as mysam
from collections import deque

def reads_alignment(ref_genome, reads_file_name, nmismatch, num_cores):
    '''
    time /exds/users/oabalbin/sw/bwa/bwa-0.5.8a/bwa aln -k 2 -t 3 
    /exds/projects/alignment_indexes/bwa/hg19_illumina/hg19_illumina.fa s_4_2_sequence.txt > 
    s_4_2_sequence.txt.sai
    '''
    outfile_name = reads_file_name.replace+'.sai'
    args=['bwa','aln' '-k', nmismatch, '-t', num_cores, ref_genome, reads_file_name]
    
    f = open(outfile_name, "w")
    retcode = subprocess.call(args, stdout=f)
    f.close
    
    return outfile_name

def paired_read_alignment_sam(ref_genome, align_reads_pair1, align_reads_pair2, raw_reads1, raw_reads2):
    '''
    time /exds/users/oabalbin/sw/bwa/bwa-0.5.8a/bwa sampe 
    /exds/projects/alignment_indexes/bwa/hg19_illumina/hg19_illumina.fa 
    ./s_4_1_sequence.txt.sai ./s_4_2_sequence.txt.sai ./s_4_1_sequence.txt 
    ./s_4_2_sequence.txt > s_4_12_sequence.aln.sam
    
    # Using the path for the gatk well formed bam files with read group, 
    and platform information
    time /exds/users/oabalbin/sw/bwa/bwa-0.5.8a/bwa sampe 
    -i MARK -m DEPRISTO -l LIBRARY -p ILLUMINA 
    /exds/projects/alignment_indexes/gatk/hg19/bwa/hg19.fa 
    ./s_3_1_sequence.hg19.sai ./s_3_2_sequence.hg19.sai 
    ./s_3_1_sequence.txt ./s_3_2_sequence.txt > s_3_12_sequence.hg19.aln2.sam
    
    '''
    fields=align_reads_pair1.split('_')
    name = fields[0]+'_'+fields[1]+'_'+'12'+fields[3]
    outfile_name=name.replace('.sai','.sam')
    args=['bwa', 'sampe', '-i', 'MARK', '-m', 'DEPRISTO', '-l', 'LIBRARY', '-p', 'ILLUMINA', 
          ref_genome, align_reads_pair1,align_reads_pair2, raw_reads1, raw_reads2]
    
    f = open(outfile_name, "w")
    retcode = subprocess.call(args, stdout=f)
    f.close
    
    return outfile_name



def main_bwa_alignment(ref_genome, read_pairs_set, bwa_param, num_cores):
    '''
    It perfomrs a paired bwa alignment
    remove duplicate, sort the reads in cordinate order, and index the sorted file
    bwa_param => dictionary with the parameters for a bea run
    nmismath= number of mismatches
    '''
    
    nmismatch = bwa_param['nmismatches']
    
    '''
    aligned_reads=deque()
    for thispair in read_pairs_file:
        alng= reads_alignment(ref_genome, thispair, nmismatch, num_cores)
        aligned_reads.append(alng)
    '''
    aligned_reads=[reads_alignment(ref_genome, thispair, nmismatch, num_cores) 
                   for thispair in read_pairs_set]

    # sampe = paired aligned file (sam), bampe paired aligned file (bam)
    # brs = bampe removed duplicates(r) and coordinate sorted (s) file
    sampe_file = paired_read_alignment_sam(ref_genome, aligned_reads[0],aligned_reads[1], 
                              read_pairs_set[0], read_pairs_set[2])
    
    '''
    bampe_file = mysam.sam2bam(ref_genome, sampe_file)
    bampe_rmdup = mysam.remove_duplicates(bampe_file)
    bampe_rmdup_sorted = mysam.sort_bam_file(bampe_rmdup)
    brs_indexed = mysam.index_bam_file(bampe_rmdup_sorted)
    '''
    # Implement that as pysam calls
    
    return brs_indexed


        
    